#include <stdio.h>
#include <math.h>
#include "winbgi2.h"

void trid(double a[], double b[], double c[], double d[], int n);
int spline(double x[], double y[], int n1, double a[], double b[], double c[], double d[]);

void main()
{
	int i, j, n, ns;
	double x[200], y[200], a[200], b[200], c[200], d[200], xx, xb, yb, xp, xk, dx,t;
//  ***** zadane wezly funkcji interpolowanej *****
           FILE *p1 = fopen("dane.txt", "r");
           fscanf(p1,"%d", &n);
    	   for (i = 1; i <= n; i = i + 1) fscanf(p1, "%lf%lf", &x[i], &y[i]);

//  ******     obliczenie wspolczynniki splinow   *******

	ns = spline(x, y, n, a, b, c, d);
	printf("wpolczynniki splinow a,b,c,d\n");
	for (i = 1; i <= ns; i = i + 1) printf("%d  %lf  %lf  %lf  %lf\n", i, a[i], b[i], c[i], d[i]);

//  *****     interpolacja,  biezace x z przedziału <xp,xk>   ****
	printf("podaj przedzial interpolacji oraz krok xp,xk,dx :  ");
	scanf("%lf%lf%lf", &xp, &xk, &dx);
	xb = xp;



}




// funkcja spline wyznacza wspolczynniki splinow   a,b,c,d
int spline(double x[], double y[], int n1, double a[], double b[], double c[], double d[])
{
	int i, n;
	double aa[200], bb[200], cc[200], dd[200];
	double ll, lp;
	//   ilosc segmentow spline'a    
	n = n1 - 1;

	//   rownanie dla brzegu lewego (warunek naturalny)

	ll = x[2] - x[1];
	bb[1] = 2 * ll;
	cc[1] = ll;
	dd[1] = 3 * (y[2] - y[1]);

	//   rownania dla punktow wewnetrznych
	for (i = 2; i <= n; i++)
	{
		ll = x[i] - x[i - 1];
		lp = x[i + 1] - x[i];
		aa[i] = lp;
		bb[i] = 2 * (ll + lp);
		cc[i] = ll;
		dd[i] = 3 * ll / lp*(y[i + 1] - y[i]) + 3 * lp / ll*(y[i] - y[i - 1]);
	}
	//   rownanie dla brzegu prawego (warunek naturalny)
	ll = x[n1] - x[n];
	aa[n1] = ll;
	bb[n1] = 2 * ll;
	dd[n1] = 3 * (y[n1] - y[n]);
	//   wyznaczenie pochodnych (nachylen), 
	//    rozwiazanie trojdiagonalnego ukladu rownan

	trid(aa, bb, cc, dd, n1);

	//   wyznaczenie wspolczynnikow spline'a
	//      y(x) = a(i)*t^3 + b(i)*t^2 + c(i)*t^1 + d(i)  ; t = x - x(i)
	for (i = 1; i <= n1; i++)
	{
		d[i] = y[i];
		c[i] = dd[i];
	}
	for (i = 1; i <= n; i++)
	{
		ll = x[i + 1] - x[i];
		b[i] = (3 / ll*(d[i + 1] - d[i]) - (c[i + 1] + 2 * c[i])) / ll;
		a[i] = (2 / ll*(d[i] - d[i + 1]) + (c[i + 1] + c[i])) / ll / ll;
	}
	return n;
}



// funkcja rozwiazywanie ukladu r-n z macierza trojdiagonalna
void trid(double a[], double b[], double c[], double d[], int n)
{
	int i;
	double q;
	for (i = 2; i <= n; i++)
	{
		q = a[i] / b[i - 1];
		b[i] = b[i] - q*c[i - 1];
		d[i] = d[i] - q*d[i - 1];
	}

	d[n] = d[n] / b[n];
	for (i = n - 1; i >= 1; i = i - 1) d[i] = (d[i] - c[i] * d[i + 1]) / b[i];
}


